module math_mod

	use const_mod

	implicit none

  private
	public solve_tridiagonal

contains

	subroutine solve_tridiagonal(coef_a, coef_b, coef_c, coef_r, x, nsize)

		!Solutions to a system of tridiagonal linear equations A*x=r
		real(r8), intent(in) :: coef_a(nsize), coef_b(nsize), coef_c(nsize), coef_r(nsize)
		real(r8), intent(out) :: x(nsize)
		integer, intent(in) :: nsize
		integer i
		real(r8) gama(nsize), rho(nsize)

		if (coef_b(1) == 0.0) then
			write(*,*) "solve_tridiagonal error"
			stop
		end if
		! step1: forward elimination
		gama(1) = coef_c(1) / coef_b(1)
		do i = 2, nsize - 1
			gama(i) = coef_c(i) / (coef_b(i) - gama(i-1) * coef_a(i))
		end do
		gama(nsize) = 0.0

		rho(1) = coef_r(1) / coef_b(1)
		do i = 2, nsize
			rho(i) = (coef_r(i) - rho(i-1) * coef_a(i)) / (coef_b(i) - gama(i-1) * coef_a(i))
		end do
		!setp2: back substitution
		x(nsize) = rho(nsize)
		do i = nsize-1, 1, -1
			x(i) = rho(i) - gama(i) * x(i+1)
		end do
		
	end subroutine solve_tridiagonal

end module math_mod